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Online  vehicle  trajectory  reshaping  is  desired  for  a  class  of  autonomous  air 
vehicles  such  as  RLVs  in  order  to  avoid  catastrophic  failure  when  subjected  to 
performance  restricting  damages  and  failures.  An  Adaptive  Trajectory  Reshaping 
and  Control 1  (ATRC)  system  is  envisioned  that  responds  to  altered  vehicle  conditions 
by  continuously  retargeting  and  reshaping  the  reference  RLV  trajectory  satisfying 
the  feasibility  constraints.  On-line  trajectory  reshaping  to  determine  a  feasible 
reference  trajectory  is  computationally  a  difficult  problem  for  real  time 
applications.  ATRC  is  exploring  the  principles  of  vehicle  dynamics  inversion  for  on¬ 
line  generation  of  feasible  reference  trajectory.  Two  essential  components  for 
generating  reference  trajectory  for  air-vehicles  using  “inverse  dynamics” 
methodology  are  aerodynamic  model  of  the  vehicle  that  is  representative  of  the 
current  state  of  the  vehicle,  and  a  framework  for  modeling  the  vehicle  trajectory. 
Physics  based  modeling  software  such  as  Missile  DATCOM  allows  fast  computation 
of  aerodynamic  coefficients  for  given  flight  points  and  the  results  can  be  stored  in 
tabular  form.  However,  for  efficient  real-time  trajectory  reshaping  application,  it  is 
desirable  to  represent  aerodynamic  coefficients  in  smooth  functional  forms  that  are 
governed  by  a  few  parameters.  Similarly,  trajectories  must  also  be  represented  by 
smooth  functions.  In  this  paper  we  present  modeling  of  smooth  functions  using  a  set 
of  basis  functions  that  are  suitable  for  aerodynamic  modeling  and  trajectory 
reshaping  of  the  air  vehicles.  A  desirable  feature  for  function  modeling  is  the  easy 
imposition  of  boundary  as  well  as  mid  point  constraints  in  the  function  using  a  small 
number  of  parameters  without  limiting  the  scope  of  the  function.  In  this  paper  we 
present  a  design  method  for  generating  orthonormal  polynomial  basis  functions  in 
one  and  two  dimensions  with  constraints. 


I.  Introduction 

The  large  potential  for  space  utilization  is  not  being  exploited  as  it  is  currently  inhibited  by  the  huge 
cost  of  launching  operations.  The  benefits  of  advanced  space  utilization  can  be  greatly  increased  by  making 
space  utilization  more  affordable  The  Reusable  Launch  Vehicle  (RLV)  programs  are  targeted  towards 
affordable  space  utilization.  However,  to  maintain  the  economical  viability  of  RLVs,  it  is  important  to 
enhance  operational  safety  and  reliability  by  providing  the  RLV  the  capability  to  respond  to  various 
uncertainties  and  evolving  emergency  situations.  Responding  to  an  uncertain  environment  after  a 
damage/failure  presents  many  tough  technical  problems  for  this  class  of  vehicles.  These  problem  manifests 


Senior  Researcher,  1408  University  Dr,  College  Station,  TX-77840,  Senior  Member  AIAA, 
averma@KBSI.com. 

fResearcher,  1408  University  Dr,  College  Station,  TX-77840,  Member  AIAA, 
kvadakkeveedu@KB  SI .  com. 

{Electronics  Engineer,  2210  Eighth  Street,  Bldg.  146,  Rm.  305,  WPAFB,  OH  45433-7531,  Member  AIAA, 
Member  AIAA,  Michael . Oppenheimer@wpafb .  af.mil 

§Senior  Aerospace  Engineer,  2210  Eighth  Street,  Bldg.  146,  Rm.  305,  WPAFB,  OH  45433-7531,  Associate 
Fellow  AIAA,  David.Doman@wpafb.af.mil 

1 

American  Institute  of  Aeronautics  and  Astronautics 


EXTENDED  ABSTRACT 


in  the  following  challenges  that  must  be  overcome:  First,  to  adequately  determine  and  model  the  dynamic 
characteristics  of  the  vehicle  in  the  altered  state  after  a  damage/failure;  second,  to  estimate  the  new 
constraints  and  limitation(s)  of  the  vehicle;  third,  to  adapt  and  reconfigure  the  command,  control  and 
guidance  of  the  vehicle  to  the  modified  system  dynamics;  and  fourth,  to  design  and  plan  a  new  feasible  path 
with  respect  to  the  end  goal  maximization.  A  high  percentage  of  such  damage/failure  cases  leave  the 
vehicle  in  an  uncontrolled  and  uncertain  environment  with  a  high  probability  of  ultimately  entering  a  state 
of  catastrophic  failure.  The  high  cost  of  loss  resulting  from  catastrophic  failures,  has  prompted  researchers 
in  the  direction  of  developing  technologies  to  assist  in  minimizing  such  failures. 

Damage  to  a  vehicle  or  a  sub-system  failure  may  result  in  modification  of  the  applicable  trajectory 
constraints  and/or  the  dynamical  behavior  of  the  system,  consequently  making  the  previously  designed 
reference  trajectory  infeasible.  An  acceptable  trajectory  for  a  dynamical  system  is  a  solution  of  a  two-point 
boundary  value  problem  for  a  set  of  governing  differential  equation  of  motion.  The  real  world  systems  such 
as  air  vehicles  are  highly  non-linear  systems  and  impose  a  set  of  constraints  on  the  trajectory  variables  as 
well  as  control  variables.  In  inverse  dynamics  approach,  a  trajectory  is  specified  first,  which  results  in 
solving  a  set  of  algebraic  equations,  yet  strictly  satisfying  the  non-linear  differential  equations  of  a  non-flat 
system. 

In  this  paper  first  we  introduce  the  architecture  of  an  Adaptive  Trajectory  Reshaping  and  Control 
(ATRC)  system  for  the  general  class  of  RLV  systems,  which  is  based  on  the  principles  as  described  in  Ref. 
[1].  Next  we  discuss  the  general  inverse  dynamics  approach  for  trajectory  reshaping  of  RLVs  and  the 
motivation  for  functional  modeling.  For  trajectory  determination,  there  are  two  types  of  functions  that  must 
be  modeled.  The  first  set  of  functions  is  needed  to  describe  the  vehicle  trajectory.  The  functions  for  the 
spatial  coordinates  of  the  vehicle  constitute  a  vehicle  trajectory.  These  trajectory  functions  are  normally 
functions  of  one  independent  variable.  The  second  set  of  functions  that  must  be  modeled  in  real  time  is  the 
modified  aerodynamic  constraints.  The  aerodynamic  constraints  are  modeled  by  defining  aerodynamic 
coefficients  for  the  flight  envelope  of  the  vehicle.  Mostly,  these  functions  are  two  dimensional  with  Mach 
and  angle  of  attack  being  the  independent  parameters.  For  faster  convergence,  it  is  desired  that  the 
functions  be  smooth,  with  continuity  in  the  value  as  well  as  first  derivative.  A  general  approach  to  model  a 
function  is  to  parameterize  the  function  and  then  determine  the  parameters  that  satisfy  any  governing 
constraints  to  a  satisfactory  level.  In  Section  IV,  we  describe  a  novel  approach  for  designing  a  set  of  basis 
functions  with  a  class  of  constraints  built  in,  that  reduces  the  complexity  of  the  solution  and  results  in  better 
approximations. 


II.  Adaptive  Trajectory  Reshaping  and  Control  (ATRC)  System 

The  ATRC  system  enhances  RLV  capability  to  avoid  catastrophic  failure  when  subjected  to 
performance  restricting  damages  and  failures.  The  overall  goal  of  ATRC  translates  into  specific 
requirements  for  design  and  development  of  functionalities  related  to  adaptable  and  reconfigurable 
command,  control,  and  guidance  system  for  the  RLVs.  and  real  time  solution  techniques  for  generating 
feasible  trajectories. 

Figure  1  shows  the  general  architecture  of  the  envisioned  ATRC  system  for  RLVs.  Note  that  the  above 
structure  is  specific  to  longitudinal  motion  of  the  vehicle;  however,  it  can  be  easily  extended  to  the  six 
degrees  of  freedom.  The  main  components  of  the  envisioned  ATRC  system  requires: 

1 .  On-line  system  identification  that  includes  physics  based  modeling  that  accounts  for  vehicle  damage, 
parameter  estimation  and  parameter  projection  for  constraint  boundary  determination.  The  constraint 
boundaries  influence  the  trajectory  reshaping  of  the  vehicle. 

2.  Real  time  trajectory  determination  for  reshaping  reference  trajectory  under  feasibility  constraints. 

3.  Adaptive,  closed  loop  control  and  guidance  system  for  reference  trajectory  tracking. 


2 

American  Institute  of  Aeronautics  and  Astronautics 


EXTENDED  ABSTRACT 


Outer  Loop  Guidance 

q  cmc. 


Tracking 

Errors 


Guidance 

Cmds. 


RLV 

Model 


ft 


Inner  Loop  Stabilization 


Adaptive 

Stabilization 


RLV 

Model 


Reference 


Trajectory 


T  raj  ectory  Re-targeting 


Reshaping 


Objs./ 

Constraints  Specs. 

* 


ATRC-Architecture 


System 


Learning  and 


Estimation 


Parameter  ID 


On 

RLV] 

[me 

Vlodd 

Model 

Learning 


Figure  1.  Architecture  of  ATRC 


III.  Inverse  Dynamics  Approach  for  Trajectory  Generation 

The  advantage  of  inverse  dynamics  approach  is  that  one  can  avoid  integration  of  differential  equations 
of  motion  altogether.  The  main  difference  between  regular  trajectory  generation  approach  and  inverse 
dynamics  approach  is  that  that  in  the  latter  approach  we  first  parameterize  the  trajectory  and  then  use 
numerical  techniques  to  solve  for  these  parameters  that  minimize  the  objective  function  and  satisfy  other 
inequality  constraints.  Note  that  once  a  smooth  trajectory  is  specified,  the  time  derivatives  of  the  trajectory 
parameters  are  also  fixed.  With  the  availability  of  state  and  its  time  derivatives,  the  differential  equations 
become  simple  algebraic  equations,  which  is  much  faster  to  solve  than  differential  equations. 

A.  Reference  Trajectory  Design 

To  determine  a  trajectory  for  a  non-linear  dynamic  system,  a  solution  must  be  found  that  satisfies  the 
set  of  differential  equations  governing  the  dynamics  of  the  system.  Further,  the  trajectory  solution  should 
not  violate  some  non-linear  constraints,  which  limits  the  operation  capability  of  the  system.  For  an  aircraft, 
the  constraints  arise  due  to  several  factors  such  as  limitations  on  the  angle  of  attack ,  load  factor , 
aerodynamic  heating ,  and  actuator  saturation. 

There  are  two  primary  approaches  for  trajectory  generation  and  these  have  been  classified  in  the 
literature  [2]  as  the  “integral  approach”  and  the  “differential  approach.”  In  any  approach,  where  generation 
of  a  trajectory  involves  the  integration  of  the  equations  of  motion,  this  approach  is  classified  as  the 
“integral  approach.”  In  a  differential  approach,  an  assumed  functional  form  for  trajectory  is  differentiated 
to  obtain  algebraic  functions  for  the  higher  derivatives,  which  are  required  to  impose  constraints  on  the 
control  inputs  for  the  “inverse  dynamics”  solution.  There  are  various  applications  where  inverse  dynamics 
have  been  used,  such  as  spacecraft  trajectories  and  path  planning  in  robotics  [3]  and  overhead  cranes  [4]. 
Historically,  the  inverse  dynamics  approach  has  been  used  for  “differentially  flat”  systems.  A  system  is 
“differentially  flat”  [see  5,6]  if  there  exists  a  set  of  outputs,  known  as  “flat  outputs,”  such  that  there  is  a 
one-to-one  correspondence  between  the  trajectories  of  flat  outputs  and  the  full  state  and  control  inputs  of 
the  system.  In  our  approach,  we  use  an  inverse  dynamics  approach  for  determining  trajectory  as  this  allows 
us  to  solve  algebraic  equations  instead  of  integrating  ODEs.  With  the  inverse  dynamics  approach  for 
aircraft  trajectories,  a  problem  arises  due  to  inherent  under- actuation  in  most  of  the  aircrafts.  For  a  six- 
degrees-of- freedom  (6DOF)  aircraft,  there  are  normally  four  controls:  thrust,  elevator,  aileron,  and  rudder. 
Ref.  [7]  defines  a  novel  trajectory  generation  scheme,  which  uses  pseudo  forces  for  inverse  dynamic 
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computation.  In  this  paper  we  concentrate  on  function  modeling  approach  to  make  the  inverse  dynamic 
solution  approach  more  efficient. 

B.  Inverse  Dynamics  Approach 

Assume  x  be  the  state  vector  for  a  vehicle  and  x  =  x(r)  represent  a  reference  trajectory  of  the  vehicle. 
The  goal  of  the  inverse  dynamic  approach  is  to  design  an  inverse  trajectory  such  that  the  trajectory  and  the 
corresponding  solution  for  control  inputs  satisfy  all  constraints.  The  governing  equations  of  motion  for  an 
air  vehicle  may  be  approximated  with  terms  that  are  non-linear  with  respect  to  some  of  the  states  and  linear 
with  respect  to  some  others.  For  example,  if  x1  and  x2  are  the  subset  of  the  state  vector  x ,  a  typical 
governing  equation  may  be  written  in  the  form 

x  =  Cj  (x1 )  +  C2  (x1  )x2  +C3(x1)m)  (1) 

where  Cx  (xx ),  C2  (xx ),  C3  (xx )  are  non-linear  functions  in  xx ,  and  the  over  all  function  is  linear  with 
respect  to  x2  and  control  inputs  u.  For  example,  the  pitch  dynamics  of  a  vehicle  for  longitudinal  dynamics 
can  be  written  as 


=  cm0  (M’ a)+  Crnq  ( M,a)q  +  Cm&  ( M,a)Se ,  (2) 

where  xl  =  \M ,ci\ ,  u  =  Se,  and  C', , C2 , C'3  are  Cm0/lYr  ,Cm  jlyy  ,Cm&  j Ijj  respectively.  The  inverse 
dynamic  solution  for  Eq.  (1)  can  be  written  as 

u  ~  Q-1  (^i  )(x  -  Cx  (xx )  +  C2  (x1  )x2 )  (3) 

However,  note  that  the  existence  of  C3”1  may  not  always  be  guaranteed  for  any  arbitrary  trajectory.  In 
the  event  of  any  damage  to  the  RLV  ,  the  reference  trajectory  must  be  re-designed  in  view  of  the  altered 
vehicle  dynamics.  For  a  feasible  solution,  a  numerical  iterative  approach  is  used  [7],  where  the  reference 
trajectory  is  perturbed  until  feasibility  is  ensured.  To  facilitate  perturbation  of  the  vehicle  trajectory,  the 
trajectory  must  be  parameterized.  Generally,  the  trajectory  design  must  meet  certain  constraints  on  end 
points,  and  sometimes  there  may  be  constraint  at  mid  points.  In  this  paper  we  will  present  functional 
modeling  that  allows  function  parameterization  as  well  as  ensures  any  constraints  on  the  reference 
trajectory.  For  fast  convergence  of  the  numerical  solution  to  the  inverse  problem  ,  it  is  desired  that  the 
vehicle  model,  represented  by  terms  Cif  z=l,--*,3be  smooth  functions.  As  the  vehicle  dynamics  are 

altered,  the  altered  coefficient  functions  Ci{xx )  must  be  re-determined  on-line. 

If  the  damage  condition  is  known,  physics  based  modeling  techniques,  such  as  the  approach  used  by 
Missile  DATCOM  [8],  can  determine  the  numerical  values  of  the  coefficient  functions  at  various  flight 
points.  In  our  approach,  we  first  determine  the  vehicle  coefficients  in  tabular  form  spanning  the  designed 
flight  envelope.  However,  we  require  a  framework  to  capture  the  coefficients  in  functional  forms  such  that 
the  non-linear  functions  are  continuous  and  smooth.  In  the  following  sections  we  address  the  design  of  a  set 
of  constrained  orthonormal  bases,  followed  by  one  dimensional  trajectory  modeling,  and  two-dimensional 
coefficient  functional  representation  in  parametric  form. 

IV.  Constrained  Orthonormal  Polynomial  Basis  Functions 

In  this  section  we  present  an  approach  to  design  orthonormal  polynomial  basis  functions  with  desired 
constraints  built  in.  The  choice  of  basis  functions  ultimately  influences  the  complexity  of  modeling  the 
desired  functions  or  behaviors.  It  is  well  known  that  a  sub  set  of  a  complete  set  of  basis  functions  can 
approximate  any  given  function  to  a  desired  accuracy  by  choosing  a  sufficient  number  of  basis  functions 
elements  in  the  sub  set.  However,  if  a  function  to  be  approximated  must  meet  certain  constraints,  it  presents 
various  problems.  First,  a  finite  set  of  basis  function  may  not  ensure  the  constraint-  satisfaction  on  the 
function  to  be  approximated.  Second,  a  large  number  of  basis  functions  may  be  required  for  the  satisfactory 
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constrained  approximation.  Now,  if  we  can  create  a  set  of  basis  functions,  where  all  the  basis  functions 
satisfy  the  given  constraints,  any  linear  combination  of  those  functions  will  also  satisfy  the  given 
constraints.  With  this  motivation  we  developed  an  approach  to  design  a  constrained  orthonormal 
polynomial  basis  functions.  To  explain  our  approach,  we  first  present  a  way  for  generating  one-dimensional 
unconstrained  orthonormal  polynomial  based  basis  functions  and  then  extend  it  to  generate  constrained 
orthonormal  polynomial  basis  functions. 

A.  One  Dimensional  Basis  Functions  with  Built-In  Constraints 

Let  ®,(x)  represent  the  zth  ID  basis  function  of  the  orthonormal  polynomial  basis  set  given  as: 

»■= ox-  (4) 

k=0,i 

Notice  that  the  zth  basis  function  has  z+ 1  polynomial  coefficients  to  be  determined.  We  impose  the  given 
constraints  on  each  of  the  basis  functions.  E.g.  if  a  desired  constraint  is  that  the  modeled  function  g(x) 
should  be  zero  at  x  =  0.5,  then  we  create  a  new  basis  function  set  where  oc  .  (x)  =  (®  (x)  (8)  (x  -  0.5)),  where 

(8)  is  the  polynomial  multiplication(or  ID  convolution  of  the  polynomial  coefficients).  g(x)  =  0  for  all  the 

functions  in  the  space  spanned  by  this  new  basis  set.  The  Gram- Schmidt  procedure  is  applied  to  this 
constrained  basis  set  to  generate  an  orthonormal  basis  set. 

The  inner  product  between  any  two  functions  in  this  space  is  defined  as 

l 

(<D (x),  ®  *  (x))  =  J  O ,  (x)®  k  (x)dx .  (5) 

0 

The  norm  and  the  orthogonality  conditions  can  be  obtained  as: 

l 

J  (®  t  ( x))2dx  =  1 ,  Normalizing  condition.  (6) 

o 

l 

J®/(x)®A:(x)A:  =  0,  k  =  Orthogonality  condition.  (7) 

o 

The  orthogonality  condition  gives  i  independent  equations  while  an  additional  equation  from  normalization 
ensures  that  all  the  coefficients  ak,k-0,i  can  be  obtained. 


For  automatically  generating  the  orthogonal  basis  set  we  use  the  Gram-Schmidt  recursive  approach. 
Assuming  that  we  already  know  <&k,k  =  0,z -1 .  Define  an  arbitrary  starting  polynomial  (x)  of  degree  z. 


For  k  =  0  to  i  -  1  we  determine  fk+x  (x)  as 


(8) 


The  zth  basis  function  <$i( x )  is  now  given  as 


(9) 


Next  we  will  extend  the  above  formulations  to  impose  a  given  set  of  constraints.  The  constraint  can  be 
specified  on  the  function  value  or  higher  derivatives  at  specified  points.  A  typical  nth  order  constraints  is 
given  as: 


dnOj(x) 

dxn 


{x=a) 


=  0 


(10) 


Before  we  describe  an  approach  to  determine  constrained  orthogonal  polynomial  basis  functions,  we 
note  some  properties  of  polynomials.  A  polynomial  satisfies  an  nth  order  constraint  as  given  in  Eq.  (10)  at 
x -a  when  a  is  the  root  of  the  reduced  polynomial  obtained  through  n  differentiations.  Also  note  that  a 
polynomial  of  degree  m  trivially  satisfy  any  constraint  of  order  more  than  m ,  while  constraints  of  order  m  or 
less  must  be  imposed  by  adjusting  the  polynomial  coefficients.  Sometimes,  a  polynomial  of  degree  m  may 
not  exist  that  satisfy  all  the  constraints,  especially  when  there  are  too  many  constraints  of  order  m  and  lees 
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that  must  be  imposed.  In  this  case  the  mth  degree  polynomial  will  not  be  a  part  of  constraint  orthogonal 
basis  functions. 

For  constrained  orthogonal  polynomial  function,  first  we  define  a  set  of  polynomial  functions  of  all 
possible  degrees  that  satisfy  the  given  constraints.  Next,  we  use  Gram-Schmidt  scheme  to  orthogonalize 
and  normalize  the  functions  to  form  a  basis  function  set.  The  ith  function  f  (x)  starting  from  i-  0  is 
determined  in  the  following  iterative  manner. 

Set  /°  (x)  =  xl . 

Next  extract  the  lowest  order  constraint  from  the  given  set  of  constraints  that  must  be  imposed.  If  the 
lowest  order  constraint  has  order  m,  and  m  <  degree  (// ) ,  then  incorporate  the  constraint.  Incorporating  a 
constraint  also  increases  the  degree  of  the  polynomial  function.  Repeat  the  process  for  the  next  lowest  order 
constraint,  until  the  lowest  order  constraint  is  equal  to  the  current  degree  of  (// ) .  Figure  2  and  Figure  3 

show  some  examples  of  a  class  of  orthogonal  basis  functions  that  incorporate  the  given  constraints.  The 
example  presents  normalized  orthogonal  basis  functions  in  the  range  0  to  1 . 


-  x  - ►  -  x  - ► 

(a)  (b) 

Figure  2.  Constrained  Polynomial  Basis  Functions,  (a)  Basis  Function  are  constrained  to  be  Zero  at 
x  =  0.5.  (b)  Both,  function  and  its  derivative  are  constrained  to  be  Zeros  atx  =  0.5 

In  Figure  2(a),  the  value  of  the  function  is  constrained  to  be  zero  mid  way,  i.e.  /(x)|x=05  =  0  Figure  2(b) 
adds  an  additional  constraint  on  the  slope  of  the  curve  or  the  first  derivative  of  the  function  given  as 
df{x)/dx\x=0  5  =  0  .  In  Figure  3,  the  application  of  the  constraint  has  been  shifted  to  the  quarter  point,  i.e. 

x  =  0.25.  In  this  case,  the  basic  nature  of  orthogonal  basis  functions  remains  the  same,  albeit  with  a  loss  of 
symmetry.  Figure  2  and  Figure  3  demonstrate  the  ability  to  generate  a  range  of  readily  available  sets  of 
constrained  orthogonal  polynomial  basis  functions. 
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-  x - ►  -  x - ► 

(a)  (b) 

Figure  3.  Constraint  Polynomial  Basis  Functions,  (a)  Basis  Function  are  Constraint  to  be  Zero  at  x  = 
0.25.  (b)  Both,  Function  and  its  Derivative  are  Constrained  to  be  Zeros  at  x  =  0.25 

B.  Two-Dimensions  Ortho-Normal  Polynomials  Basis  Functions 

In  this  section  we  demonstrate  an  approach  to  produce  a  set  of  polynomial  orthonormal  basis  functions 
in  two  dimensions,  represented  by  X  and  y  with  maximum  degree  of  N ,M  respectively.  First  we 

introduce  the  vector  space  of  two-dimensional  polynomials  and  then  the  algorithm  to  generate  orthonormal 
basis  vectors. 


f(x  =  .25)  =  0, 
f\x  =  .25)  =  0 


7.  Matrix  Representation  of  2D  Polynomials 

Let  PN>M  be  a  set  of  polynomials  where  P(x,y)s  PN>M  is  a  polynomial  in  x  and  y  of  degree  N  and  M 
respectively.  Let  the  coefficients  be  aif . 


P(x,y)=aNMxN yM  +aN_mxN  xyM  +---  +  a0Ux°yM  +aN_XMxN  lyM  +•••  +  %//, 


or 


cV-i  .  yM~J 

Let  us  represent  the  polynomial  coefficients  in  matrix  form  as, 


aNM 

aN-\M 

aN-2M 

a0  M 

A  = 

aNM-\ 

aN-2M-\ 

a0  M-] 

l 

o 

aN- 10 

aN- 2,0 

o 

o 

^3 

JM+l,N+l 


Then  the  polynomial  P{x,y )  can  be  represented  in  matrix  form  as, 

P(x,y)  =  XT AY , 

where. 


(11) 


(12) 


(13) 
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~xN  ~ 

V 

x  = 

xN~x 

,  and 

7  = 

yM-l 

1 

N+\ 

l 

(14) 


2.  Representation  of  a  2D  Polynomial  Inner  Product  Vector  Space 

Let  PN>M  be  a  set  of  polynomials  in  x  and  y  of  degrees  up  to  N  and  M  respectively.  The  polynomial  set 
P N)M  defines  a  polynomial  vector  space  over  R  (set  of  real  numbers)  with 

®  =  K;  I  ®i,j  =  X  yJ  lo<(<7V,0< j<M 

as  the  basis  vectors  that  satisfy  the  basic  properties  of  associativity,  commutativity,  existence  of  identity 
and  inverse,  and  properties  of  scalar  multiplication.  In  order  to  define  orthogonality,  we  must  define  an 
inner  product  space.  An  inner  product  space  is  a  vector  space  with  an  additional  structure  of  inner  product. 
We  define  the  inner  product  on  the  2D  polynomial  vector  space  as, 

1  i 

( P(x,y),Q(x,y))  =  ^P(x,y)®Q(x,y)dxdy ,  P(x,y),Q(x,y)  e PNM.  (16) 

oo 

Our  goal  is  to  construct  an  orthonormal  basis  for  the  polynomial  vector  space  PN>M  .  As  a  first  step  we 

compute  an  orthonormal  basis  functions  set  of  vectors  using  the  Gram-Schmidt  procedure.  Then  we 
approximate  the  function  based  on  its  projection  on  the  basis  functions. 

3.  Constrained  Orthonormal  Basis  Set  for  2D  Functions 

Many  times  we  would  like  to  represent  a  class  of  functions  in  terms  of  basis  functions  that  must  satisfy 
certain  constraints.  For  an  efficient  representation,  we  would  like  to  construct  a  set  of  basis  functions 
satisfying  the  given  constraints  in  addition  to  the  orthonormality  condition.  In  this  report  we  restrict  our 
discussion  to  three  types  of  constraints;  constraints  on  (i)  the  function  value ,  (ii)  function  value  and  first 
derivative  or  (iii)  only  on  the  first  derivative.  In  the  following  we  describe  the  generation  of  2D 
constrained  basis  functions. 

(i)  If  the  constraints  are  imposed  only  on  the  value  of  the  function  along  the  curve  C(x,y)  =  0,  the 
polynomial  basis  set  0  is  modified  to  introduce  the  constraints  and  create  a  new  basis  set  0C . 
The  new  basis  set  is  given  by, 

©c  =  {0ij{x,y)®c{x,  y)\. 

The  symbol  ®  denotes  polynomial  multiplication  (or  equivalently  discrete  convolution  of  the 
polynomial  coefficients).  Any  function  g(x,y)  that  lies  in  the  space  spanned  by  0C  satisfy  the 
constraint  that,  along  the  curve  C(x,y)  =  0,  g(x,y)  =  0. 

Note  that  the  basis  set  0C  is  not  yet  orthonormal.  The  constrained  orthonormal  basis  set  Oc  is 
constructed  using  the  Gram-Schmidt  orthogonalization  procedure.  See  Figure  4  for  an  example  of 
the  orthonormal  basis  functions  generated  using  this  procedure.  The  basis  set  0C  and  the 
orthonormal  basis  set  ®c  span  the  same  function  space. 
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(ii)  If  the  constraints  are  imposed  on  the  value  as  well  as  the  derivatives  along  a  curve  C(x,y)  =  0, 
a  constrained  basis  set  0C  is  computed  from  the  original  basis  set  0  as, 

0c  =  ^ij{x,y)®c(x,y)®c{x,y)). 

The  symbol  (x)  denotes  polynomial  multiplication  (or  equivalently  discrete  convolution  of  the 
polynomial  coefficients).  Any  function  g(x,y)  that  lies  in  the  space  spanned  by  0C  satisfy  the 

constraint  that,  along  the  curve  C(x,y)  =  0,  g(x,y)  =  0  and  g'(x,y)  =  0.  The  constrained  orthonormal 
basis  set  ®c  is  constructed  using  the  Gram-Schmidt  procedure  as  discussed  in  the  previous 
section.  Figure  5  shows  an  example  of  the  orthonormal  basis  functions  generated. 

(iii)  Constraint  only  on  the  derivative  of  the  function.  We  then  modify  the  given  polynomial  0  as, 

©c  =  KoE  {©-  KoEAl®  c(x,y)®c(x,y) . 

Figure  6  gives  an  example  of  orthonormal  functions  generated  with  constraints  only  on  the  derivatives. 


Figure  4.  Orthonormal  basis  function  set  of  polynomials  in  x  and  y  of  degrees  2  and  2,  with  the 
constraint  that  at  x+y  =  1,  the  value  of  the  function  is  zero. 
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Figure  5.  Orthonormal  basis  function  set  of  polynomials  in  x  and  y  of  degrees  2  and  2,  with  the 
constraint  that  at  x+y  =  1,  the  value  and  derivatives  of  the  function  are  zero. 


Figure  6.  Orthonormal  basis  function  set  of  polynomials  in  x  and  y  of  degrees  2  and  2,  with  the 
constraint  that  at  x+y  =  1,  the  derivative  of  the  function  is  zero. 
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V.  One-Dimensional  Function  Modeling  with  Constraints  using  Basis  Functions 

Often  we  need  to  estimate  or  model  non-linear  functions  that  must  satisfy  some  constraints.  For 
example,  the  inverse  dynamic  solution  approach  requires  parameterization  of  multiple  system  output 
variables  or  the  trajectory  variables  with  some  constraints  imposed  on  them.  Note  that  the  feasible 
trajectory  must  also  satisfy  the  governing  differential  equations.  The  trajectory  parameters  are  determined 
using  an  optimization  process  that  minimizes  the  error  between  trajectory  variables  and  the  solution  of  the 
governing  equations.  Note  that  if  some  of  the  trajectory  constraints,  such  as  the  boundary  constraints  are 
forced  in  the  trajectory  parameterization  itself,  it  helps  in  reducing  the  complexity  of  the  optimization 
problem.  This  work  extends  the  approach  followed  in  Verma  et  al  [7,9]  for  representing  trajectory  in 
functional  form  that  is  useful  to  impose  boundary  and  in-point  constraints.  To  prescribe  a  smooth  trajectory 
functional  for  a  specific  trajectory  coordinate,  the  position  coordinate  is  represented  by  a  twice 
differentiable,  smooth  function  so  that  velocity  and  acceleration  can  be  uniquely  defined.  The  path  for  an 
individual  position  coordinate  is  defined  as  a  function  of  normalized  time  t  -t  IT f  ,  where  Tf  is  the  total 

time  for  the  maneuver.  An  individual  coordinate  trajectory  is  structured  to  have  two  parts,  a  base  trajectory 
and  a  perturbation  of  the  base  trajectory.  If  P{r )  is  the  ith  coordinate,  it  is  chosen  to  be  of  the  form 

P(r)=P0(r)+Pl(r)^«/^/(r)  (17) 

j= i 

Here  Pq(t)  is  the  base  trajectory  function,  which  is  chosen  such  that  it  satisfies  any  boundary  or  mid 
point  constraints  for  the  trajectory.  One  example  approach  P0(r)  could  be  chosen  as  a  minimum  degree 

polynomial  spline  that  satisfies  the  desired  constraint  conditions.  The  base  trajectory  by  itself  may  not  be  a 
true  solution  of  the  governing  equations  and  hence  an  infeasible  trajectory.  To  make  the  trajectory  feasible, 
we  add  a  second  term  on  the  right  hand  side  of  the  trajectory  equation.  The  perturbation  term  must  be 
designed  to  ensure  that  the  overall  trajectory  function  P{y)  is  the  solution  of  the  governing  equations. 
Perturbation  term  is  defined  using  a  set  of  orthonormal  basis  functions  \cp.  (r)\,j  -  1,-  •  •,/  that  are  modified 

by  the  term  Px  (r) .  The  term  Px  (r)  is  a  weight  polynomial  constraining  the  perturbation  term  to  contribute 
zero  value  for  the  already  satisfied  desired  constraints  by  base  trajectory  P0(r).  For  example  Px{r)  can  be 
chosen  as 

- a)p {z -b)q  ,  (18) 

where  p  and  q  are  integers,  usually  less  or  equal  to  3.  The  integers  p  and  q  depend  upon  the  highest 
order  condition  match  required  at  times  r  =  a  and  r  =  b  respectively.  For  example,  when  boundary 
conditions  at  z  =  0  and  r  =  1  up  to  acceleration  level  must  be  imposed  on  the  trajectory  function,  Px{r)  can 

be  chosen  as  r3(r-l)3 .  The  drawback  of  this  approach  is  that  the  after  basis  functions  (pj  are  modified  by 

Px  (r) ,  it  looses  the  orthogonality,  which  makes  it  computationally  harder  to  determine  the  coefficients  for 
the  trajectory  solution. 

In  this  paper,  we  extend  the  one  dimensional  function  modeling  by  incorporating  the  constraints  in  the  basis 
functions  set  as  explained  in  the  Section  IV.A,  thus  eliminating  the  need  for  enforcing  constraints  by  an 
extra  term  Px  (r) . 


VI.  Aerodynamic  Coefficient  Function  Estimation 

An  important  goal  of  the  ATRC  system  is  the  adaptive  reshaping  of  the  RLV  trajectory  in  the  presence 
of  altered  dynamic  characteristics  of  the  vehicle  when  unexpected  damage  occurs  in  the  various  operating 
scenarios.  Any  damage  to  a  vehicle  that  has  an  impact  on  the  external  shape  of  the  vehicle,  or  that  creates 
an  impediment  in  normal  functioning  of  the  control  surfaces,  results  in  alteration  of  the  vehicle’s 
aerodynamic  characteristics.  Figure  7  and  Figure  8  show  few  examples  of  the  Pitch  moment  coefficient 
variation  in  the  presence  of  various  damage  scenarios.  Since  the  aerodynamic  behavior  of  the  vehicle  is 
captured  in  aerodynamic  coefficients  that  are  used  for  the  design  of  vehicle  control  and  trajectory  planning, 
it  becomes  mission  critical  to  adapt  reference  trajectory  for  the  altered  vehicle  dynamics.  Hence  we  need  an 
approach  to  build  a  smooth  on-line  aerodynamic  model.  Physics  based  modeling  software  such  as  Missile 
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DATCOM  allows  fast  computation  of  aerodynamic  coefficients  for  given  flight  points  and  the  results  can 
be  stored  in  tabular  form.  However,  for  efficient  real-time  trajectory  reshaping  application,  it  is  desired  to 
represent  aerodynamic  coefficients  in  smooth  functional  form  that  are  governed  by  few  parameters.  In  this 
section,  we  present  a  piecewise  continuous  and  smooth  function  model  in  two  dimensions. 


Figure  7.  Moment  Coefficient  with  AOA  for 
Nominal  and  Various  Failure  Cases 


Figure  8.  Moment  Coefficient  with  Mach  for 
Nominal  and  Various  Failure  Cases 


A.  Physics  Based  Aerodynamic  Modeling 

Given  the  geometry  of  the  vehicle,  very  good  estimations  of  the  aerodynamic  coefficients  can  be 
generated  based  on  Physics  based  modeling  in  a  very  small  time.  For  example,  if  damage  on  the  vehicle  is 
characterized  adequately,  the  DATCOM  technology  allows  specification  of  the  altered  geometry  of  the 
vehicle  to  compute  the  corresponding  aerodynamic  coefficients  at  desired  points  of  the  flight  envelope.  In 
our  approach,  we  generate  a  large  number  of  data  points  for  the  aerodynamic  coefficients  using  DATCOM 
technology  and  then  fit  a  piecewise  continuous  and  smooth  function  to  create  a  functional  model  for 
vehicle  aerodynamics.  The  functional  form  is  later  used  for  trajectory  reshaping. 

B.  Finite  Element  Function  Approximation 

We  formulate  the  aerodynamic  coefficient  function  with  a  set  of  parameters.  First  we  show  our 
approach  for  piecewise  continuous  and  smooth  one-dimensional  function,  which  is  later  extended  to  two- 
dimensions.  We  used  a  finite  element  modeling  approach  so  that  the  approximation  function  would  capture 
local  variations  in  an  efficient  manner.  Ref.  [10,  11]  demonstrates  the  use  of  finite  element  piecewise 
approximation  for  mapping  geopotential.  Ref.  [6,  8]  applied  the  technique  for  aerodynamic  coefficients 
representation.  First,  the  argument  space  of  the  function  is  divided  to  form  a  grid  with  one  control  point  at 
every  grid  point.  At  each  control  point  we  use  a  local  polynomial  function  that  is  determined  using  a 
weighted,  least  square  method  from  a  given  set  of  nominal  data  generated  from  DATCOM. 

If  a  function  is  available  at  few  data  points,  a  continuous  function  can  be  obtained  by  interpolation. 
Note  that  interpolation  in  input  space  of  9lp  is  a  surface  in  ${p+l  space.  The  interpolation  in  strict  sense 
can  be  defined  in  the  following  manner.  Suppose,  given  N  points:  qi,q2,mm '  ,q^  P  and  N  scalars: 
ax ,  a2 ,  •  •  • ,  aN  e  91 ,  we  wish  to  find  a  function  /  such  that  / (qi)-ai,  i  =  f  •  •  • ,  N .  A  standard  way  of 

interpolation  is  to  use  the  Lagrange  polynomials.  However,  when  number  of  data  points  is  large,  the  degree 
of  polynomial  becomes  high  that  may  result  in  excessive  undesirable  modulations  in  the  functions.  Another 
simple  approach  is  to  linearly  interpolate  between  the  adjacent  points.  The  linear  interpolation  preserve  the 
function  values  at  the  given  data  points,  but  smoothness  is  lost.  For  smoothness,  the  continuity  of  the  first 
derivative  of  the  approximated  function  must  be  ensured.  In  some  cases,  not  only  there  is  more  knowledge 
about  the  function  value,  but  also  we  may  know  or  we  can  estimate  the  slope  of  function  behavior.  For 
example,  the  DATCOM  model  can  estimate  the  aerodynamic  coefficients  and  certain  dynamic  derivatives 
at  a  given  flight  point.  In  this  case  it  is  desired  that  besides  continuous  derivative  for  smoothness,  the  values 


12 

American  Institute  of  Aeronautics  and  Astronautics 


EXTENDED  ABSTRACT 


of  first  derivatives  are  also  preserved.  In  our  functional  modeling  approach  for  vehicle  aerodynamic  model, 
we  wish  to  preserve  the  function  value  and  the  estimated  function  slopes. 

4.  One-Dimensional  Finite  Element  Approximation 

A  major  benefit  of  using  a  finite-element  approach  is  that  in  this  approach  local  approximations  for  the 
function  model  are  used  that  are  simpler  and  more  accurate.  Here  we  present  a  finite  element  approach  for 
modeling  the  function.  First  we  define  a  set  of  control  points  that  divides  the  input  domain  in  many  finite 
elements.  Next,  at  each  control  point  we  obtain  a  local  polynomial  approximation  of  the  function  based 
upon  the  input  data.  The  scope  of  each  local  polynomial  centered  at  the  control  point  lies  in  between  the 
adjacent  grid  points  (see  Figure  9).  Notice  the  overlapping  of  local  polynomials  that  helps  in  obtaining  a 
smooth  function  over  the  entire  range.  Once  local  approximations  are  determined,  a  smooth  global 
approximation  is  obtained  as  a  weighted  combination  of  the  local  approximations.  The  smoothness  of  the 
approximation  function  implies  that  the  function,  as  well  as  its  first  derivative,  is  continuous.  Notice  that 
grids  are  not  restricted  to  be  equidistant.  To  capture  nonlinearity  effectively,  more  control  points  should  be 
placed  near  high  non-linearity.  Figure  9  uses  a  second  order  polynomial  function  for  local  approximations. 

A  highlight  of  this  Finite  Element  approximation  is  that  it  preserves  the  local  function  value  and  its  first 
derivative  at  its  control  point.  This  is  achieved  by  a  smooth  weighting  function  that  gracefully  goes  from 
unity  to  zero,  from  one  control  point  to  another,  without  contributing  to  a  first  derivative  at  both  control 
points.  The  weighting  function  for  the  two  overlapping  local  approximation  curves  is  given  as 

Wx  (q)  =  (1  -  q  f  (1  +  2 q) ,  W2(q)  =  \-  Wx (q)  , 

where  q  is  the  normalized  coordinate  with  the  origin  at  the  first  control  point.  Note  that  the  weighting 
function  ensures  that  Wl(0)  =  l,  W2( 0)  =  0  at  the  first  control  point  and  similarly,  Wl(l)  =  0 ,  W2(l)  =  l 
at  the  second  control  point.  If  fx(q )  and  f2(c j)  are  two  local  approximations,  the  weighted  global 
approximation  F(q)  between  two  control  points  is  given  as 

F(q)  =  Wx  (q)fi  (q)  +  W2  (q)f2  (q) . 

Figure  10  shows  the  final  approximated  function  that  is  made  of  the  weighted  combinations  of  the  local 
function  approximations. 


Figure  9:  Local  Approximations  Centered  at 
Control  Points 


Figure  10:  Smooth  Functional  Approximation 
using  Weighted  Local  Approximations 


5.  Multi-Dimensional  Finite  Element  Approximation 
In  this  section  we  extend  the  formulation  for  FEM  modeling  for  multi  dimensional  functions.  Let  F  be 
given  as 
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F  =  F(q1,q2,---,q„) ,  (19) 

and  we  have  to  determine  an  estimate  F(qx,q2,-  •  •  ,qn)  from  a  finite  element  model.  As  in  one¬ 
dimensional  case,  we  assume  that  the  domain  of  F  is  covered  with  finite  number  of  equidistant  nodes. 
Each  preliminary  local  approximation  if  .  ... .  is  valid  in  the  2x2x---x2  hypercube  centered  at  a  node 

represented  by  n  indices  given  as  (z1?/2?‘“?4)  •  Ea°h  index  is  numbered  in  increasing  order  for  the 

nodes  lying  in  its  dimension.  Consider  a  hyper  cube  formed  by  2”  nodes,  each  node  representing  a  comer 
of  the  hyper  cube.  The  final  approximation  if  .  .../  is  valid  in  the  unit  hypercube  whose  “lower  left 

corner”  is  (/j ,  i2 ,  •  •  • ,  in  )  .  We  shift  the  origin  of  the  coordinates  to  the  node  (z^ ,  i2 ,  •  •  • ,  in  )  that  has  the 
lowest  indices  in  each  dimension.  The  local  coordinates  are  normalized  as 


9j3iJ+i-qUj 

Let  us  define  a  function  FT  as  a  function  of  the  normalized  coordinate  gv  as 

W  (Fj )  =  (1  -  #/  )2  (1  +  2(jj ) 

Now  the  weight  function  for  an  arbitrary  comer  node  (z^ ,  i2 ,  •  •  • ,  ik  + 1,  ik+l  + 1, 


(20) 

(21) 

in  )  of  the  hyper  cube 


can  be  written  as 

"U.-A+  '  ^(1  -  ^)^(1  -  ^+1)  • 

? 

for  0<qj  <1,  j  =  b--n  ^ 

The  weight  function  defined  by  Eq.  (22)  has  the  following  properties: 

•  Weight  is  unity  at  its  own  node  and,  it  is  zero  at  any  other  node  or  edge  formed  by  other 
nodes. 

•  Summation  of  all  the  weights  corresponding  to  the  comer  nodes  of  the  hyper  cube  at  any 
point  inside  the  hyper  cube  is  equal  to  one. 

The  estimated  function  F  can  be  defined  in  terms  of  normalized  coordinates  centered  at  node 


(  h  5  4  ?  *  *  *  ?  4  ) 


as 


F(qx,q2 


w 


j=l  .  (23) 

6.  An  Example  of  Finite  Element  Function  Approximation  for  Functions  of  Two  Variables. 

Let  Z  be  the  matrix  of  values  of  a  function  sampled  at  the  grid  points  of  a  2D  space  defined  by  x  andy.  We 
are  interested  in  representing  the  function  Z  in  closed  form,  using  2D  polynomial  approximation.  The 
closed  form  representation  can  be  used  to  evaluate  the  function  Z  at  arbitrary  x  and  y  points.  We 
approximate  the  function  Z  using  a  number  of  2D  polynomial  function  elements  of  a  smaller  area  of 
support.  The  number  of  function  elements  to  be  used  and  the  degrees  of  the  polynomials  in  x  and  y  are 
inputs  to  the  algorithm.  We  specify  this  by  the  x  andy  coordinates  of  the  control  points  as  separate  vectors. 
Note  that  the  control  points  lie  exactly  on  the  grid  points.  At  these  data  points,  the  approximated  value  of 
the  function  Z  equals  the  value  of  the  local  function  centered  at  that  control  point.  The  chosen 
approximation  is  locally  optimal  in  the  least  squared  error  criteria  (i.e.,  within  the  area  of  support  of  each 
function  element). 

To  formulate  finite  element  function  approximation  satisfying  certain  constraints,  we  proceed  as  follows. 
The  orthonormal  basis  set  of  functions  are  computed  that  satisfy  the  constraints  and  then  searched  for  local 
functions  that  are  optimal  in  the  vector  space  spanned  by  the  constrained  orthonormal  basis  functions, 
based  on  least  square  optimality  criteria.  More  details  on  the  finite  element  function  approximation  can  be 
found  in  [7]. 


Concentrating  on  the  finite  element  function  approximation  method  for  2D  functions,  we  initially 
proceeded  by  calculating  the  weights  in  the  x  and  y  directions  separately.  The  actual  weight  function  at 
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each  point  in  the  grid  is  the  product  of  the  x  and  y  weights  at  that  location.  The  x  and  y  weight  functions 
satisfy  the  criteria  that  at  the  current  control  point  location,  it  has  a  value  of  one  and  gradually  decreases  to 
zero  as  it  moves  away.  The  products  of  the  ID  weight  functions  also  satisfy  the  same  criteria  in  the  2D  x-y 
plane.  The  intended  consequence  of  this  property  is  that  the  approximated  function  value  will  be  equal  to 
the  value  of  the  local  function  at  the  control  point. 

Figure  1 1  shows  an  example  of  a  weight  function  in  the  2D  plane  that  is  the  product  of  the  two  ID  weight 
functions. 


Weight  function  centered  at  a  control  point  (i  j) 


Figure  11.  Weight  function  centered  at  a  control  point.  The  value  is  one  at  the  location  of  the  control 
point  and  gradually  decreases  and  reaches  zero  at  the  neighboring  control  points. 

Figure  12  shows  the  results  of  finite  element  approximation  for  the  CM  coefficient  for  a  non-damaged 
(baseline)  spacecraft  with  16  control  points  along  the  Mach  axis  and  21  control  points  along  the  alpha  axis. 
The  basis  functions  were  polynomials  of  degree  2  in  x  and  2  in  y.  . 


Given  function  -  Cm  :  Baseline  FEM  approximation  *  Cm-  Baseline 


FEM  Error  Baseline 


Figure  12.  Finite  element  function  approximation,  (a)  Given  function  the  CM  coefficient  for  baseline 
(b)  FEM  Approximated  Function  (c)  Error  between  the  given  and  approximated  functions.  16 
control  points  along  the  Mach  axis  and  21  control  points  along  the  alpha  axis  were  used  for  FEM 
approximation.  Each  local  function  is  a  polynomial  of  degree  2  in  x  and  2  in  y. 
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VII.  Conclusion 

In  this  paper  we  addressed  the  function  modeling  techniques  to  facilitate  an  efficient  on-line  inverse 
dynamics  methodology  for  trajectory  reshaping  of  the  RLVs.  We  presented  a  technique  for  generating 
orthonormal  polynomial  basis  functions  with  in-built  function  constraints  at  divers  points  of  the  function. 
The  general  approach  allows  the  constraints  ranging  from  the  function  value,  to  constraints  on  values  of  the 
higher  derivatives  of  the  functions.  Next  we  presented  a  finite  element  approach  for  modeling 
multidimensional,  piecewise  continuous  and  smooth  functions  The  use  of  FEM  based  smooth  modeling  and 
use  of  constrained  orthonormal  basis  function  are  expected  to  improve  the  convergence  of  iterative 
numerical  solution  for  dynamic  inversion  problem. 
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